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ABSTRACT :  Reliable  predictions  of  relative  binding  free  energies  are  essential  in  drug  discovery,  where  chemists  modify  promising 
compounds  with  the  aim  of  increasing  binding  affinity.  Conventional  thermodynamic  integration  (TI)  approaches  can  estimate 
corresponding  changes  in  binding  free  energies  but  suffer  from  inadequate  sampling  due  to  the  ruggedness  of  the  molecular  energy 
surfaces.  Here,  we  present  an  improved  TI  strategy  for  computing  relative  binding  free  energies  of  congeneric  ligands.  This  strategy 
employs  a  specific,  unphysical  single-reference  (SR)  state  and  Hamiltonian  replica  exchange  (HREX)  to  locally  enhance  sampling. 
We  then  apply  this  strategy  to  compute  relative  binding  free  energies  of  12  ligands  in  the  L99A  mutant  of  T4  lysozyme.  Besides  the 
ligands,  our  approach  enhances  hindered  rotations  of  the  important  VI 11  as  well  as  V87  and  LI  18  side  chains.  Concurrently,  we 
devise  practical  strategies  to  monitor  and  improve  HREX-SRTI  efficiency.  Overall,  the  HREX-SRTI  results  agree  well  (R2  =  0.76, 
RMSE  =  0.3  kcal/mol)  with  available  experimental  data.  When  optimized  for  efficiency,  the  HREX-SRTI  precision  matches  that  of 
experimental  measurements. 


1.  INTRODUCTION 

High-quality  predictions  of  binding  free  energies  of  small 
molecules  to  their  biomolecular  targets  are  important  in  drug 
design.  The  continued  growth  of  computational  power  has 
enabled  applications  of  statistical  mechanics-based  free  energy 
perturbation  (FEP)  and  thermodynamic  integration  (TI) 
methods  to  real  life  problems.1-'^  These  advanced  computa¬ 
tional  techniques  are  considered  gold  standards  for  binding 
free  energy  predictions,  akin  to  isothermal  titration  calorimetry 
(ITC) — a  technique  that  measures  the  binding  free  energies 
experimentally.6,7 

Since  their  inception,  great  progress  has  been  made  in  improv¬ 
ing  FEP  and  TI  methods.  Introduction  of  the  soft-core  potentials 
has  made  the  calculations  more  reliable.8  10  Subsequently, 
multiconfiguration  simulation  protocols11  have  laid  the  founda¬ 
tion  for  the  application  of  generalized  ensemble  strategies  such 
as  Hamiltonian  replica  exchange  (HREX)12  17  which  further 
improve  the  quality  of  the  FEP  and  TI  simulations.18-22  Mean¬ 
while,  better  postprocessing  protocols  have  been  developed 
that  resulted  in  more  reliable  predictions  of  free  energies  and 
assessments  of  corresponding  standard  deviations.4, 23- 

Despite  this  progress,  binding  free  energy  calculations  remain 
challenging  because  of  sampling  limitations  that  are  inherent  in 
the  molecular  dynamics  (MD)  methods  used  in  the  simulations. 
Conventional  FEP  and  TI  free  energy  calculations  are  known  to 
be  sensitive  to  starting  conformations  of  the  bound  complexes. 
For  example,  in  a  well-studied  L99A  mutant  of  T4  lysozyme,  the 
conformation  of  the  binding  site  residue  VI 11  affects  binding 
free  energies  of  indene  and  p-xylene  ligands  by  as  much  as  6  kcal/ 
mol.1,28-30  Other  hindered  residues  can  have  a  similar  effect.1 
The  predictions  also  depend  on  the  initial  orientation  of  the 
ligands  in  the  binding  pocket.  These  challenges  are  likely  to  be 
general  and,  therefore,  need  to  be  properly  addressed. 
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Many  enhanced  sampling  approaches  have  been  devised  to 
combat  conformational  challenges.  It  is  impossible  to  list  all  of 
them  here,  but  we  will  name  a  few  that  benefit  alchemical  free 
energy  calculations.  One  of  the  earliest  approaches  scaled  parts  of 
the  potential  energy  before  and  after  an  alchemical  transforma¬ 
tion  to  enhance  sampling.31  Recently,  an  approach  called  accel¬ 
erated  molecular  dynamics  (AMD)  was  combined  with 
alchemical  free  energy  calculations.32-34  This  approach  adds  a 
boosting  potential  to  reduce  barriers  and  is  also  independent  of 
alchemical  transformations.3'^  Other  methods  exploited  the 
alchemical  transformations  to  enhance  sampling.  Examples  in¬ 
clude  A  dynamics36  38  or  its  Monte  Carlo  (MC)  counterparts, 
such  as  chemical  MC/MD39’40  and  more  general  simulated 
scaling.41,42  In  fact,  these  methods  share  important  features  with 
currently  developing  FEP  and  TI  methods  augmented  with 
Hamiltonian  replica  exchange  (HREX). 15-22,43 

Enhanced  sampling  will  overcome  dependence  of  FEP  and 
TI  predictions  on  the  starting  conformations  and  ultimately 
improve  their  accuracy  and  precision.  Judiciously  combining 
multiconfigurational  FEP  and  TI  methods  with  HREX  can  achieve 
enhanced  sampling.2-’43  Previously,  we  presented  a  TI  variant, 
called  single-reference  TI  augmented  with  HREX  (HREX-SRTI), 
which  achieved  convergence  of  solvation  free  energies  for  a 
challenging  case  of  an  amide  system  where  conventional  FEP 
and  TI  methods  failed.22  The  amide  in  question  had  an  internal 
barrier  to  cis/trans  interconversion  that  was  insurmountable  in 
conventional  MD  simulations.  HREX-SRTI  was  able  to  generate 
converged  results  using  simulation  times  of  only  4  ns. 

In  the  present  study,  we  apply  HREX-SRTI  to  the  well-studied 
T4  lysozyme  mutant.  Although  this  system  has  a  simple  binding 
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site — a  hydrophobic  cavity  buried  beneath  the  protein  surface — 
it  is  sufficiently  complex  to  render  conventional  FEP  and 
TI  approaches  ineffective.  Importantly,  our  binding  free 
energy  predictions  for  this  system  can  be  compared  to  the  pre¬ 
viously  published  independent  computational  and  experimental 
values.1’  8-30 

First,  we  assess  the  variability  in  the  free  energy  predictions 
using  regular  SRTI.  Then,  we  employ  SRTI  with  the  HREX 
option  and  demonstrate  that  HREX  efficiency  is  crucial  to 
obtaining  converged  results.  Thus,  we  provide  practical  recipes 
to  improve  the  efficiency  of  HREX-SRTI  simulations.  Finally, 
we  use  one  of  these  recipes  to  optimize  HREX  efficiency  and 
obtain  highly  converged  results  for  the  most  challenging  of  the 
ligands — indole. 


2.  METHODS 

2.1.  Parameters  for  Small  Molecule  Ligands.  We  stud¬ 
ied  12  small  molecules  as  follows:  benzene,  toluene,  o-xylene, 
p-xylene,  ethyl-benzene,  n-propyl-benzene,  n-butyl-benzene,  ;-butyl- 
benzene,  phenol,  indene,  benzofuran,  and  indole.  The  initial  coordi¬ 
nates  of  all  of  the  small  molecules  in  the  present  study  were  derived 
using  the  program  CORINA44  Where  applicable,  we  performed 
conformational  expansion  using  the  program  ROTATE.45  Two 
charge  models  described  in  sections  2.1.1  and  2.1.2  were  used  in 
combination  with  the  GAFF  force  field  for  small  molecules.46  Each 
conformational  ensemble  was  structurally  refined  through  geometry 
optimization  using  the  AMI  semiempirical  quantum  mechanical 
potential47  as  implemented  in  MOPAC7,  version  1.1 1.4S 

2.1.1.  AM1BCC  Charge  Model.  Partial  charges  from  each 
unique  conformation  were  accumulated  using  Boltzmann  weight¬ 
ing  with  appropriate  degeneracies  by  their  AMI  energies  at  300  K. 
The  final  AMI  charges  were  symmetrized  where  applicable  and 
then  augmented  through  the  BCC  procedure47,50  implemented  in 
the  ANTECHAMBER  program51'52  from  AMBER  TOOLS,  ver¬ 
sion  1.2.  The  resulting  conformation-independent,  properly  sym¬ 
metrized  set  of  AM1BCC  charges  is  expected  to  reproduce  HF/ 
6-31G(d)  RESP  charges  to  a  good  approximation.5'  56 

2.1.2.  RESP  Charge  Model.  For  a  select  subset  of  molecules, 
we  derived  HF/6-31G(d)  restricted  electrostatic  potential  fit 
(RESP)  partial  charges  using  B3LYP/6-31G(d)  optimized  geom¬ 
etries.  The  geometry  optimization  and  the  electrostatic  potential 
calculations  were  achieved  using  Gaussian  09.57  The  RESP  fit58 
was  performed  using  the  ANTECHAMBER  program. 

2.2.  Setup  of  Protein-Ligand  Complexes.  Coordinates  of  all 
of  the  ligands  in  complexes  with  the  L99A  mutant  of  T4  lysozyme 
were  derived  from  a  crystal  structure  (PDB:  181L)  of  the  protein 
bound  to  benzene.57  Initial  placements  of  ligands  other  than  benzene 
were  derived  using  graph  theory.  Specifically,  molecules  were  repre¬ 
sented  as  graphs  on  the  basis  of  their  atom  and  bond  types. 
Subsequently,  association  graphs  were  constructed,  and  maximal 
cliques  were  found  to  match  atoms  in  the  benzene  rings  of  each 
molecule  to  those  of  the  bound  benzene.60  There  are  12  different 
cliques  that  give  rise  to  12  unique  placements  of  each  ligand  within  the 
binding  site  of  the  protein.  Some  of  the  cliques  are  degenerate 
depending  on  the  symmetry  of  the  molecule.  Thus,  for  benzene,  all  1 2 
cliques  are  degenerate,  yielding  identical  complex  structures  that  differ 
only  in  the  numbering  of  the  carbon  atoms  of  the  ligand.  However,  for 
ligands  such  as  benzofuran,  indene,  and  indole,  each  clique  yields  a 
unique  conformation  of  the  complex 

By  construction,  initial  protein  coordinates  and  those  of  the 
benzene  ring  are  identical  for  all  of  the  complexes,  while  the 
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Figure  1.  Thermodynamic  cycle  for  computing  binding  free  energies 
relative  to  an  unhindered,  unphysical  reference  state  using  SRTI.  Hor¬ 
izontal  arrows  represent  alchemical  transformations  of  indole  into  the 
benzene  core  in  water  and  the  protein  binding  site.  The  protein  residue 
VI 11  is  alchemically  converted  to  p-A  at  the  same  time.  In  the  reference 
state,  the  molecular  volumes  of  the  ligand  and  the  VI 11  residue  (shown 
with  green  and  brown  shaded  contours,  respectively)  are  reduced  due  to 
disappearing  atoms.  The  disappearing  atoms  interact  with  the  rest  of  the 
system  through  soft-core  Coulomb  and  vdW  potentials  and  are  connected 
to  the  reference  core  by  dashed  bonds.  The  torsional  potentials  associated 
with  the  virtual  atoms  are  removed.  These  changes  greatly  enhance 
translational,  rotational,  and  torsional  degrees  of  freedom  involved  in 
the  alchemical  transformation  with  the  help  of  HREX,  thus  activating 
rotation  of  the  VI 11  and  of  the  ligand  in  the  binding  site. 


atoms  protruding  from  the  benzene  ring  change  their  position.  In 
cases  with  branched  ligands,  such  as  n-propyl-,  n-butyl-,  and 
i-propyl-benzene,  the  protrusions  have  been  examined  for  steric 
clashes  with  protein  side  chains. 

The  protein  is  described  by  an  all-atom  Amber  99SB  molecular 
mechanics  force  field  compatible  with  GAFF.61  The  solvation 
effects  were  modeled  using  a  cubic  periodic  box  of  explicit  TIP3P 
water  molecules  that  extended  at  least  10  A  beyond  the  solute.  The 
protein  system  was  neutralized  by  adding  nine  Cl-  ions. 

2.3.  Single  Reference  State.  The  choice  ofthe  reference  states 
in  SRTI  determines  which  degrees  of  freedom  of  the  system  would 
be  accelerated.22  To  allow  for  enhanced  sampling  of  the  ligands  in 
the  confines  of  the  binding  pocket,  we  chose  the  benzene  core 
without  hydrogen  atoms  as  the  reference  state  (Figure  l).  Con¬ 
veniently,  our  ligand  reference  state  is  independent  of  the  charge 
model  because  the  benzene  core  atoms  have  no  charges  by 
construction.  Previously,  we  successfully  employed  this  ligand 
reference  state  in  computing  free  energies  of  hydration.22 

Similarly,  choosing  an  appropriate  protein  reference  state 
could  enhance  sampling  of  the  hindered  protein  side  chains  such 
as  Val,  Leu,  lie,  and  Thr.  Each  hindered  residue  could  be  mutated 
to  a  modified  Ala  residue  referred  to  as  pseudo-Ala  (p-Ala). 
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Table  1.  The  L99A  T4  Lysozyme  Mutant  Relative  Binding  Free  Energies  (Standard  Deviations)  to  a  Series  of  Benzene 
Derivatives" 


compound 

exptl 

FEP 

SRTI 

Diff, 

HREX-SRTI 

Diff2 

Wat% 

Prt% 

AM1BCC  Charge  Model 

benzene 

0.0(0.2) 

0.0 

0.0(0.3) 

0.0 

0.0(0.2) 

0.0 

34 

29 

phenol 

2.5(N/A) 

3.3 

2. 2(0.6) 

-0.3 

23(0.4) 

-0.2 

25 

26 

toluene 

-0.3(0.2) 

0.0 

—0. 1(0.3) 

0.2 

-0.2(0.2) 

0.1 

30 

27 

ethylbenzene 

— 0.6(0.2) 

-1.8 

— 0.2(0.5) 

0.3 

— 0.9(0.6) 

-0.3 

27 

25 

n-propylbenzene 

-1.4(0.2) 

-1.3 

—0.4(0. 5) 

0.9 

—  1. 5(0.4) 

-0.2 

26 

23 

n-butylbenzene 

—  1. 5(0.2) 

-0.3 

—  1. 0(1.2) 

0.5 

— 1.4(1. l) 

0.1 

23 

21 

i-butylbenzene 

— 1.3(0.2) 

-0.5 

—  1.4(1. 2) 

-0.1 

-13(0.7) 

0.1 

24 

21 

o-xylene 

0.6(0.2) 

3.3 

1. 0(0.4) 

0.4 

0.8(03) 

0.2 

28 

25 

p-xylene 

0.5  (0.2) 

1.0 

0.5  (0.6) 

0.0 

0.4(03) 

-0.1 

29 

26 

indene 

0.1  (0.2) 

2.8 

1.3(0.4) 

1.3 

1.4(03) 

1.3 

21 

20 

indole 

0.3(0.2) 

4.1 

2. 9(1. 7) 

2.6 

2. 8(0.5) 

2.5 

14 

21 

benzofuran 

— 0.3(0.2) 

1.0 

0.4(0. 5) 

0.6 

0 .2(0.4) 

0.5 

24 

22 

RMS  (cyclic) 

(1.0) 

1.7 

(0.4) 

1.7 

RMS  (acyclic) 

(0.7) 

0.5 

(0.5) 

0.2 

RMS  (all) 

(0.8) 

1.0 

(0.5) 

0.9 

RESP  Charge 

Model 

indene 

—0.1  (0.4) 

-0.2 

0.0(03) 

-0.1 

b 

20 

indole 

0.6(1. 9) 

0.3 

11(1-8) 

0.8 

b 

19 

benzofuran 

— 0.4(0.6) 

-0.1 

-0.8(0.4) 

-0.5 

b 

22 

RMS  (cyclic) 

(0.8) 

0.2 

(0.7) 

0.5 

RMS  (all)c 

(1.1) 

0.4 

(1.1) 

0.3 

11  Energies  are  in  kcal/mol  relative  to  benzene.  Averages  and  standard  deviations  are  over  eight  (cyclic)  or  four  (acyclic)  independent  simulations  with 
distinct  starting  positions  of  the  ligand  in  the  binding  site  of  the  protein  and  two  (all)  independent  simulations  in  water.  Experimental  and  previously 
reported  absolute  binding  free  energies  for  benzene  are  —5.2  and  —4.6  kcal/mol,  respectively.  Because  the  effect  of  HREXon  the  simulations  in  water  is 
small,  regular  simulations  were  performed.  Diff,  is  the  free  energy  difference  between  regular  SRTI  and  experimental  results;  Diff2  is  the  difference 
between  HREX-SRTI  and  experimental  results.  The  Wat  and  Prt  columns  contain  acceptance  ratios  for  the  corresponding  legs  of  the  thermodynamic 
cycle  (Figure  l).  For  these  calculations,  p  =  2,  and  default  values  of  (X  =  1.5  parameters  were  used.  Each  simulation  was  run  using  12  TI  windows  equally 
spaced  in  X  and  4-ns-long  MD  runs  in  the  NPT  ensemble  at  1  atm  and  300  K.  The  terms  exptl  and  FEP  represent  previously  published  experimental  and 
computational  free  energy  perturbation  benchmarks.30  c  Combined  values  with  RESP(3>HF/6-31G(d)  for  cyclic  and  AM1BCC  for  acyclic  compounds. 


In  p-Ala,  hydrogen  atoms  of  the  methyl  side  chain  are  united  with 
the  C/3  carbon.  A  mutation  of  a  hindered  residue  to  p-Ala  in  the 
reference  state  would  render  its  side  chain  atoms  starting  with  C/3 
virtual.  Thus,  HREX  would  enhance  hindered  rotations  about 
bonds  such  as  Cot— C/3,  C/3— Cy,  and  outward.  In  addition,  the 
alchemical  mutation  of  binding  site  residues  to  p-Ala  can  render 
the  binding  site  in  the  reference  state  bigger,  further  aiding  in  the 
sampling  of  ligand  transitions. 

In  this  study,  we  simultaneously  enhance  sampling  of  the 
ligand  and  the  protein  side  chains  by  combining  the  correspond¬ 
ing  unphysical  reference  states. 

2.4.  TI  Simulation  Setup.  In  order  to  run  the  simulations, 
we  employed  GROMACS  version  4.0.5  in  single  precision. 
Because  the  simulated  system  is  described  by  Amber  99SB61 
and  GAFF46  molecular  mechanical  force  fields  that  are  not  native 
to  GROMACS,  we  used  the  PERL  conversion  script,  which  was 
described  previously.62  The  script  also  automates  setup  of  the 
alchemical  transformation  from  the  real  to  the  reference  state. 

2.4.1.  Soft-Core  Potentials.  In  order  to  avoid  the  end-point 
catastrophe  at  the  reference  state,  we  employed  soft-core  elec¬ 
trostatic  and  LJ  potentials9,10  as  implemented  in  GROMACS.63-67 
Earlier  calculations  employed  a  GROMOS  style  soft-core  po¬ 
tential  (eqs  1  —  3).  Here,  A  is  the  Hamiltonian  coupling  para¬ 
meter,  p  is  the  coupling  power,  r  is  the  distance  between  a  given 


pair  of  atoms,  a  is  the  soft-core  parameter,  and  O  is  the  radius  of 
interaction  computed  from  LJ  parameters.  For  certain  polar 
hydrogen  atoms,  O  is  undefined,  and  in  those  cases  a  fixed  value 
is  used.  Originally,  we  used  p  =  2,  a  =  1.5,  and  O  =  0.3  as 
recommended  in  the  user  manual.  Subsequently,  to  improve  the 
acceptance  ratio  and  level  its  distribution  over  TI  window  pairs, 
we  used  an  alternative  soft-core  potential  with  p  =  l.10  For  the 
latter  potential,  we  reoptimized  the  value  of  a  to  arrive  at  a  =  0.4. 
The  optimization  of  the  (X  parameter  was  performed  to  achieve 
the  best  convergence  behavior  by  monitoring  standard  devia¬ 
tions  across  independent  TI  runs.  Other  more  complicated 
measures  could  be  used  to  search  for  better  alchemical  paths, 
but  these  were  not  pursued  in  this  study.68 

V£W  =  (1  ~X)  VA(RA(r,X))  +  XVB(RB(r,X))  (l) 
RA(r,X)  =  (cloaXp  +  r6)1/6  (2) 

Rs(r,  X)  =  (aa6B(l-A)p  +  r6)1/6  (3) 

2.4.2.  MD  Simulation  Parameters.  The  production  runs  were 
performed  in  the  NPT  ensemble  at  T  =  300  K  and  P  =  1  atm, 
following  the  equilibration  protocol  described  previously.22 
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The  production  ran  employed  a  Langevin  thermostat  and  a 
Berendsen  barostat,64  67  with  identical  collision  frequencies  of  2  ps  1 . 

Throughout  the  simulations,  all  of  the  bonds  containing 
hydrogen  atoms  were  constrained  using  LINCS,69  and  the 
integration  time  step  was  set  to  2  fs.  We  employed  the  particle 
mesh  Ewald  (PME)  approach  for  electrostatics64-67  with  a  1  nm 
real  space  cutoff  and  switched  off  van  der  Waals  interactions  over 
the  range  of  0.8— 0.9  nm.  Typically,  production  runs  were  2-ns- 
long  for  each  TI  window,  but  in  some  cases  they  were  extended 
to  4  ns.  The  coordinates  of  the  system  were  recorded  every  1000 
steps  for  subsequent  analyses. 

2.4.3.  Regular  SRTI  Simulations.  To  obtain  the  alchemical  free 
energies  or  reversible  works,  the  real  and  reference  states  of  each 
system  corresponded  to  values  0  and  1  of  the  Hamiltonian 
coupling  parameter  X,  respectively.  Each  alchemical  SRTI  trans¬ 
formation  employed  M  equally  separated  X  windows.  The 
majority  of  the  simulations  used  M  =  12,  but  in  some  cases 
simulations  were  performed  with  M  =  23.  All  X  windows  of  an 
SRTI  simulation  had  the  same  initial  configuration.  For  each  X 
window,  we  recorded  dV/dX  values  at  every  time  step.  The  mean 
values  (dV/dX)  for  all  of  the  X  windows  were  assembled  into  the 
final  work  using  the  Fourier  beads  integration  procedure,  which 
was  described  previously.22 

Averages  of  the  final  work  values  and  their  standard  deviations 
were  computed  using  several  independent  simulations.  Specifi¬ 
cally,  for  proteins  in  complex  with  cyclic  molecules  (benzofuran, 
indene,  and  indole),  we  performed  eight  simulations  each  with 
distinct  starting  positions  of  the  ligand.  For  acyclic  (not  cyclic) 
molecules,  we  only  performed  four  such  simulations.  Finally,  for 
all  of  the  ligands  in  water,  we  performed  two  independent 
simulations.  The  differences  between  the  alchemical  work  values 
in  water  and  protein  environments  yielded  the  relative  binding 
free  energies  with  respect  to  the  unphysical  reference  state.  The 
final  relative  binding  free  energies  and  their  standard  deviations 
were  reported  with  respect  to  benzene. 

2.4.4.  HREX-SRTI  Simulations.  In  order  to  ran  HREX-SRTI 
simulations,  we  employed  an  in-house  PERL  script  interfaced 
with  GROMACS.  Replica  exchanges  were  attempted  every  1000 
MD  steps  or  2  ps.  For  the  majority  of  the  simulations,  we 
attempted  exchanges  a  total  of  2000  times  resulting  in  4-ns-long 
simulations  of  each  window.  In  special  cases,  the  number  of 
exchange  attempts  was  reduced  to  1000,  decreasing  the  simula¬ 
tion  time  to  2  ns  per  window.  Following  the  exchanges,  each  X 
window  received  a  new  random  seed  to  restart  its  MD  ran.  All  of 
the  other  simulation  details  were  the  same  as  those  associated 
with  regular  SRTI  simulations. 

2.4.5.  Analysis  of  SRTI  and  HREX-SRTI  Results.  The  analysis  of 
real  state  trajectories  was  performed  with  standard  GROMACS 
tools.  Specifically,  the  g_angle  program  was  used  to  obtain  time 
series  of  dihedral  angles  of  the  hindered  side  chains.  For  the  Val 
side  chain,  we  gathered  data  on  the  Ha— Ca—  C/3— H/3  (k1) 
dihedral,  and  for  the  Leu  side  chain,  data  on  Ha— Ca— C/3— C y 
(k  | )  and  Ca—  C/3— Cy— Hy  ( k2 )  dihedrals  were  collected.  For 
HREX-SRTI  simulations,  the  real  state  trajectory  had  to  be 
assembled  from  short  trajectories  using  an  in-house  PERL  script 
that  followed  the  state  through  all  of  the  exchanges. 

3.  RESULTS  AND  DISCUSSION 

In  order  to  demonstrate  the  utility  of  the  SRTI  approach  in 
computing  relative  binding  free  energies,  we  studied  ligand  bind¬ 
ing  to  a  well-defined  binding  site  in  the  L99A  mutant  of  T4 


lysozyme.  Specifically,  we  chose  a  congeneric  series  of  12  ligands 
derived  from  benzene  that  has  been  studied  previously.28-30 
We  distinguished  two  classes  of  compounds  within  the  series 
according  to  their  structure  outside  the  common  benzene  motif, 
namely  cyclic  and  acyclic.  Thus,  benzofuran,  indene,  and  indole 
were  considered  cyclic,  whereas  all  of  the  remaining  compounds 
were  considered  acyclic.  Because  all  of  the  ligands  in  the  series 
(Table  l),  with  the  exception  of  benzene,  can  have  multiple 
orientations  in  the  binding  pocket,  this  system  presents  a  con¬ 
siderable  challenge  for  binding  free  energy  calculations. 

The  cyclic  and  acyclic  ligands  behave  differently  when  in 
complex  with  the  protein.  For  acyclic  ligands,  the  benzene  ring 
can  flip  without  an  overall  structural  change  to  the  complex.  For 
cyclic  ligands,  the  benzene  flip  alters  the  overall  structure  of  the 
complex.  Hence,  for  cyclic  ligands,  i.e.,  benzofuran,  indene,  and 
indole,  we  selected  eight  orientations  (four  for  each  of  the  two 
states  resulting  from  the  benzene  flip).  The  degeneracy  with 
respect  to  the  benzene  flip  allowed  us  to  reduce  the  number  of 
representative  orientations  to  four  for  the  remaining  acyclic 
ligands.  Thus,  for  each  ligand,  we  performed  simulations  with 
different  initial  orientations  in  the  binding  pocket. 

The  conformations  of  the  active  site  residues  are  equally 
important  and  should  be  considered  when  determining  binding 
free  energies.28  30  For  the  L99A  mutant  of  T4  lysozyme,  the 
conformational  state  of  the  VI 11  side  chain  profoundly  affects 
the  computed  binding  free  energies.  Because  this  residue  lines  the 
surface  of  the  binding  site,  inadequate  sampling  of  its  conforma¬ 
tions  has  been  shown  to  cause  discrepancies  of  as  much  as  6  kcal/ 
mol.28-30  Other  residues  in  the  active  site  may  have  similar  effects 
on  binding  free  energies.1  Hence,  we  need  to  improve  sampling  of 
the  ligand  and  relevant  protein  conformations  at  the  same  time. 
This  combination  makes  the  problem  particularly  challenging. 

3.1 .  Regular  SRTI  Simulations.  3.1.1.  AM  1 BCC Charge  Model. 
The  average  variability  in  the  relative  binding  free  energies  from 
regular  SRTI  is  0.7  kcal/mol  (Table  l).  Most  of  the  ligands  have 
standard  deviations  in  the  range  of  0.3— 0.6  kcal/mol.  However, 
the  largest  acyclic  ligands,  n-butyl-benzene  and  i-butyl-benzene, 
show  increased  standard  deviations  of  1.2  kcal/ mol.  Surprisingly, 
one  of  the  cyclic  ligands,  indole,  exhibits  a  record  high  standard 
deviation  of  1.7  kcal/mol.  Indole,  like  other  cyclic  ligands,  is 
expected  to  have  hindered  flip  transitions  in  the  binding  site. 

In  comparison  to  experimental  values,  the  regular  SRTI 
approach  has  an  RMSE  of  1.0  kcal/mol  relative  to  and  excluding 
benzene.  Smaller  ligands  are  in  good  agreement — within  0.3 
kcal/mol — with  the  available  experimental  data  (Table  l).  It 
should  be  noted  that  only  an  upper  estimate  of  binding  free 
energy  is  available  for  phenol,  which  does  not  bind  the  T4 
lysozyme  mutant  well.  For  more  extended  acyclic  molecules, 
such  as  n-propyl-  and  n-butyl-benzenes,  the  agreement  is  not  as 
favorable,  with  n-propyl-benzene  demonstrating  the  highest 
deviation  of  0.9  kcal/mol.  The  predicted  binding  free  energy 
for  i-butyl-benzene  is  serendipitously  within  0.1  kcal/mol  of  the 
experimental  value.  Interestingly,  the  cyclic  ligands  exhibit  the 
largest  deviations  of  all  compounds,  diverging  by  as  much  as  2.6 
kcal/ mol  in  the  case  of  indole. 

The  disagreement  between  computed  and  experimental  bind¬ 
ing  free  energies  for  the  heterocyclic  compounds  is  instructive.  In 
particular,  the  results  for  indene  and  benzofuran  ligands  were 
well-converged  judging  by  the  low  standard  deviations  (Table  l). 
In  contrast,  the  binding  free  energy  of  indole  exhibited  a  large 
standard  deviation  of  1.7  kcal/mol.  These  observations  suggest 
that  issues  other  than  sampling  could  be  responsible  for  the 
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Figure  2.  Side  chain  conformational  transitions  for  the  indole  complex 
with  the  L99A  mutant  of  T4  lysozyme  using  the  V87p-A:Vlllp-A: 
LI  18p-A  reference.  Time  series  and  respective  histograms  are  presented 
for  the  K  |  (Ha—  Ca— C/3— H/3)  torsion  of  VI 11  and  V87  and  for  the  K1 
(Ha— Ca— C/3— Cy)  and  k2  (Ca-C/3-Cy-Hy)  torsions  ofL118. 
The  bins  of  the  histograms  were  5°  wide.  The  time  series  reports  the 
corresponding  values  of  the  torsions.  The  top  panel  summarizes  the 
regular  SRTI  simulation  results  with  relatively  few  conformational 
transitions,  whereas  the  bottom  panel  shows  the  HREX-SRTI  results 
with  numerous  such  transitions.  The  reported  HREX-SRTI  simulations 
employed  12  windows  and  the  optimized  p  =  1  soft  cores. 


Figure  3.  Side  chain  conformational  transitions  for  indole  complex  with 
the  L99A  mutant  of  T4  lysozyme  using  the  Vlllp-A  reference.  Time 
series  and  respective  histograms  are  presented  for  the  Kl 
(Ha—  Ca— C/5— H/3)  torsion  of  VI 1 1.  The  bins  of  the  histograms  were 
5°  wide.  The  time  series  report  the  corresponding  values  of  the  torsions. 
The  top  panel  summarizes  the  regular  SRTI  simulation  results  with 
relatively  few  conformational  transitions,  whereas  the  bottom  panel 
displays  the  HREX-SRTI  results  with  numerous  such  transitions.  The 
reported  HREX-SRTI  simulations  employed  12  windows  and  the 
optimized  p  =  1  soft  cores. 


overall  disagreement  with  experimental  results.  Because  we 
employed  the  AM1BCC  charge  model,  which  approximates 
the  RESP  HF/6-31G(d)  partial  charges,  we  decided  to  assess 
the  effect  of  the  model. 

3. 1.2.  RESP  Charge  Model  for  Cyclic  Compounds.  Using  RESP 
HF/6-31G(d)  charges  for  the  cyclic  compounds  considerably 
improved  the  agreement  of  regular  SRTI  predictions  with  experi¬ 
mental  values.  Indeed,  while  simulations  with  the  AM1BCC 
charge  model  systematically  overestimated  the  binding  free 
energies,  with  the  RESP  charges  the  disagreement  was  no  longer 
systematic  and  remained  within  0.3  kcal/ mol.  Thus,  for  indole, 
the  RESP  charges  lowered  the  disagreement  with  experimental 
values  by  more  than  2  kcal/mol.  Despite  the  improved  accuracy 
of  the  predictions,  the  RESP  charges  did  not  affect  convergence 
of  the  cyclic  ligands.  Indole  still  had  the  largest  standard  deviation 
of  1.9  kcal/mol.  These  results  suggest  that  a  charge  model  can 
strongly  affect  the  accuracy,  but  not  necessarily  the  precision,  of 
the  computed  binding  free  energies. 

Large  standard  deviations  in  computed  binding  free  energies 
identify  ligands  most  sensitive  to  the  initial  complex  config¬ 
uration.  The  reasons  for  this  sensitivity  likely  reside  in  hin¬ 
dered  ligand  motions.  Indeed,  inspecting  trajectories  of  the 
ground  state  simulations  for  benzofuran  using  all  12  starting 
configurations,  we  found  that  they  converge  to  only  four 


metastable  configurations.  These  metastable  configurations  de¬ 
monstrate  that  the  cyclic  ligands  do  not  flip  freely  in  the  binding 
pocket  on  a  nanosecond  time  scale.  Furthermore,  alignment  of 
the  ligands  suggests  that  the  binding  pocket  has  the  geometry  of  a 
flattened  prolate  ellipsoid.  Each  unique  configuration  contributes 
distinctly  to  the  binding  free  energy,  thus,  increasing  the  standard 
deviation.  Therefore,  enhanced  sampling  of  the  ligands  in  their 
complexes  would  allow  the  metastable  configurations  to  rapidly 
interconvert,  ultimately  improving  the  quality  of  predictions.  Because 
the  regular  SRTI  approach  does  not  have  the  ability  to  enhance  these 
transitions  alone,  we  need  to  invoke  the  HREX  option. 

3.2.  Improving  the  SRTI  Results  with  HREX.  The  use  of 
SRTI  with  HREX“2  could  simultaneously  enhance  sampling  of 
the  ligand  and  select  protein  side  chains.  First,  we  attempted 
to  enhance  motions  of  the  ligand  and  the  VI 11  side  chain. 
Specifically,  by  choosing  the  benzene  core  as  a  ligand  refer¬ 
ence  state,  we  intended  to  activate  rotations  of  the  bulkier 
ligands  within  the  active  site.  Furthermore,  by  choosing  a 
VI 1  lp-A  mutant  as  a  protein  reference  state  (see  the  Methods 
section  for  definition),  we  expected  to  activate  the  VI 11 
rotations  that  have  activation  barriers  in  the  5  —  8  kcal/mol 
range.29,70  This  strategy  is  expected  to  reduce  the  size  of  each 
ligand  to  the  size  of  the  benzene  core  while  simultaneously 
enlarging  the  binding  site  around  p-Alll  (Figure  l)  to 
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Table  2.  Improving  HREX-SRTI  Efficiency  for  the  L99A  T4 
Lysozyme  Mutant  in  Complex  with  Indole" 


system 

Nsrti 

p 

a 

AGprt(SD),  kcal/mol 

NhrEX 

accept, 

V87p-A:Vl  1  lp-A:L118p-A 

SRTI 

12 

2 

1.5 

34.2(3.8) 

SRTI 

23 

2 

1.5 

36.9(3.3) 

HREX-SRTI 

23 

2 

1.5 

39.5(1.4) 

1000b 

30 

HREX-SRTI 

12 

1 

0.3 

29.4(0.9) 

2000 

17 

HREX-SRTI 

12 

1 

0.4 

30.6(0.8) 

2000 

13 

HREX-SRTI 

12 

1 

0.5 

33.6(1.2) 

2000 

9 

Vlllp-A 


SRTI 

12 

2 

1.5 

—0.6(l.8) 

HREX-SRTI 

12 

2 

1.5 

—  1.0(l.8) 

2000 

19 

HREX-SRTI 

23 

2 

1.5 

—  1. 6(0.5) 

2000 

48 

HREX-SRTI 

12 

1 

0.4 

-1.7(0.2) 

2000 

29 

“  These  results  are  representative  of  the  protein  leg  of  the  thermody¬ 
namic  cycle  (Figure  l)  using  RESP(3>HF/6-31G(d)  point  charges  on  the 
ligand.  WSRTI  and  NHrex  refer  to  the  number  of  TI  windows  and 
exchange  cycles,  respectively.  Soft-core  parameters  involved  in  optimi¬ 
zation  are  p  and  CL  (see  text  for  description).  Averages  and  standard 
deviations  (SD)  are  over  eight  independent  simulations  with  distinct 
starting  positions  of  the  ligand  in  the  binding  site.  Two  unphysical 
references,  namely,  V87p-A:Vlllp-A:L118p-A  and  Vlllp-A,  are  con¬ 
sidered.  By  design,  the  use  of  the  former  reference  with  HREX  should 
enhance  side  chain  torsions  of  the  V87,  VI 11,  and  LI  18  along  with 
rotation  and  flipping  of  the  ligand.  The  latter  reference  should  enhance 
torsions  of  the  VI 11  side  chain  and  rotations  of  the  ligand.  Unless 
otherwise  stated,  all  TI  windows  were  run  in  the  NPT  ensemble  at  1  atm 
and  300  K  for  4  ns.  b  Each  window  was  run  for  2  ns. 


activate  ligand  flip  transitions  which  are  important  in  the 
evaluation  of  cyclic  compounds. 

Interestingly,  regular  SRTI  simulations  suggest  that  the  VI 1 1 
rotation  does  not  significantly  affect  relative  free  energies  of 
many  ligands  in  the  series.  The  free  energy  barrier  for  the  valine 
rotation  is  such  that  it  could  spontaneously  rotate  on  a  time  scale 
of  several  nanoseconds.22'70  Indeed,  we  observed  that  during 
4-ns-long  unenhanced  SRTI  simulations  of  the  ground  state, 
VI 11  does  spontaneously  flip  a  few  times  without  having  a 
significant  effect  on  the  binding  free  energies  of  many  molecules 
in  the  series  (Figures  2  and  3,  top  panels). 

The  absence  of  the  previously  reported  effect  of  VI  n  28-30  in 
cases  of  p-xylene  and  indene,  among  other  ligands,  is  not 
surprising.  It  could  be  explained  by  the  fact  that  SRTI  computes 
relative  binding  free  energies  as  opposed  to  the  absolute  free 
energies  reported  previously.  Indeed,  in  the  present  SRTI  setup, 
we  do  not  need  to  sample  an  empty  binding  site  of  the  protein, 
where  VI 11  has  a  preferred  conformation.  9  Nevertheless,  for 
the  largest  acyclic  ligands,  the  effect  could  still  be  significant. 
Therefore,  we  expect  that  enhancing  the  rotation  of  the  VI 11 
side  chain  through  HREX  would  improve  the  overall  agreement 
with  experimental  data.  This  approach  can  also  be  applied  to 
enhance  VI 1 1  rotation  in  the  absolute  binding  free  energy  calcula¬ 
tions  with  TI. 

3.2. 1.  AM1BCC  Charge  Model  and  HREX-SRTI.  Simply  turning 
on  the  HREX  option  in  SRTI  with  the  VI 1  lp-A  reference  state 
produced  seemingly  modest  improvements  over  the  correspond¬ 
ing  regular  SRTI  results  (Table  l).  The  average  standard 
deviation  for  all  of  the  ligands  decreased  from  0.7  to  0.5  kcal/ 
mol  with  the  AM1BCC  charge  model.  However,  some  of  the 


ligands  enjoyed  significantly  lower  standard  deviations  versus 
regular  SRTI.  For  example,  i-butyl-benzene  and  indole,  which 
had  among  the  highest  standard  deviations,  experienced  signifi¬ 
cant  drops  from  1.2  to  0.7  kcal/mol  and  from  1.7  to  0.5  kcal/ mol, 
respectively. 

For  the  AM1BCC  model,  the  overall  agreement  with  experi¬ 
mental  values  improved  only  slightly,  with  RMSE  decreasing  to 
0.9  kcal/mol.  Most  of  the  improvement  was  achieved  for  the 
acyclic  ligands,  which,  when  separated  from  the  rest  of  the 
ligands,  showed  a  change  in  RMSE  from  0.5  to  0.2  kcal/mol. 
The  cyclic  compounds  with  AM1BCC  charges  are  unaffected  by 
HREX  and  persistently  show  an  RMSE  of  1.7  kcal/mol.  Despite 
the  improved  agreement  with  experimental  data,  larger  acyclic 
ligands  still  exhibit  elevated  standard  deviations,  particularly  in 
the  case  of  M-butyl-benzene. 

The  lack  of  significant  improvements  could  indicate  that  other 
residues  in  the  binding  site  might  be  important.  Most  notably,  all 
of  these  simulations  used  12  windows  and  consequently  had 
modest  acceptance  ratios  (Table  l).  Indole  in  water  had  the 
lowest  acceptance  ratio  of  14%,  which  increased  to  21%  in  the 
protein  environment  with  the  Vlllp-A  reference.  Indole  and 
phenol  are  the  only  two  ligands  in  the  series  that  have  polar 
hydrogen  atoms  and  an  increased  acceptance  ratio  in  the  protein 
environment. 

3.2.2.  RESP  Charge  Model  and  HREX-SRTI.  Different  charge 
models  behave  distinctly  when  running  SRTI  with  the  HREX 
option.  Thus,  HREX-SRTI  with  RESP  charges  for  the  cyclic 
compounds  diminished  the  agreement  with  experimental  results 
compared  with  that  in  regular  SRTI.  Specifically,  RMSE  for  the 
three  cyclic  compounds  increased  from  0.2  to  0.5  kcal/mol. 
However,  the  average  standard  deviations  decreased  from  0.9  to 
0.8  kcal/mol  (Table  l).  Unexpectedly,  indole  exhibited  a  sharply 
increased  standard  deviation  of  1.8  kcal/mol  using  the  RESP 
model  compared  to  0.5  kcal/mol  with  the  AM1BCC  model. 

3.3.  Improving  HREX-SRTI  Predictions.  To  understand  the 
reason  for  poor  convergence  of  the  free  energy  values  with 
HREX-SRTI  and  to  possibly  improve  the  results,  we  examined 
the  indole  system  in  greater  detail.  Specifically,  we  focused  on  the 
indole  simulations  with  RESP  charges  that  exhibited  the  largest 
disagreement  with  experimental  values  and  the  largest  standard 
deviation. 

3.3. 1.  Triple  Mutant  Reference  State  V87p-A:Vl  I  lp-A:L  1 18p- 
A.  We  hypothesized  that  even  with  the  HREX  option  turned  on, 
flipping  of  the  indole  might  still  be  impeded  in  the  Vlllp-A 
reference  state.  In  order  to  test  this  hypothesis,  we  created  a  triple 
mutant  reference  state  by  mutating  two  additional  residues,  V87 
and  LI  18  to  p-A.  These  two  residues  pin  the  benzene  moiety  of 
the  ligands  to  the  floor  of  the  binding  site.  It  should  be  noted  that 
LI  18  has  been  experimentally  shown  to  exhibit  conformational 
variability  similarly  to  VI 11. 1 

The  triple  mutant  system  V87p-A:Vlllp-A:L118p-A  should 
completely  remove  the  ligand  flipping  restriction.  Moreover,  this 
unphysical  reference  state  could  potentially  open  water  access  to 
the  hydrophobic  binding  site  of  the  L99A  T4  lysozyme.  In  the 
proposed  experiment,  three  residues  would  undergo  alchemical 
transformations,  making  30  atoms  of  the  protein  (nine  for  each 
valine  and  12  for  leucine)  disappear  in  the  reference  state.  While 
this  should  theoretically  enhance  sampling  of  the  three  residues 
along  with  the  ligand,  it  might  be  difficult  in  practice  to  achieve  a 
sufficient  overall  exchange  rate  in  HREX-SRTI.22  Hence,  this  test 
would  also  identify  the  limits  of  our  approach  in  extending  it  to 
concurrent  activation  of  multiple  residues. 
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The  triple  mutant  V87p-A:Vl  1  lp-A:Ll  18p-A  presents  a  chal¬ 
lenge  for  HREX-SRTI.  As  seen  in  Table  2,  regular  SRTI 
simulations  with  12  and  23  windows  had  standard  deviations 
of  3.8  and  3.3  kcal/mol,  respectively,  for  just  the  protein  leg  of  the 
thermodynamic  cycle  (Figure  l).  The  fact  that  the  standard 
deviation  with  the  triple  mutant  reference  is  more  than  twice  as 
large  as  that  with  the  single  mutant  one  (1.8  kcal/mol)  suggests 
that  conformations  of  residues  V87  and  LI  18  may  indeed  affect 
ligand  binding.  The  standard  deviation  was  reduced  to  1.4  kcal/ 
mol  when  employing  23  windows  in  the  HREX-SRTI  simula¬ 
tions.  The  corresponding  12-window  simulations  did  not  achieve 
a  sufficient  overall  exchange  rate  to  yield  results  that  were  distinct 
from  the  regular  SRTI  simulations. 

3.3.2.  Assessing  HREX  Efficiency.  3.3.2. 1.  Round-Trip  Count. 
Achieving  efficient  exchanges  is  critical  for  obtaining  converged 
free  energies  in  HREX  simulations.43  The  overall  acceptance 
ratio  might  be  inadequate  to  assess  the  efficiency  of  HREX 
simulations.  While  the  overall  acceptance  ratio  is  a  satisfactory 
criterion  to  assess  the  efficiency  of  the  more  widely  applied 
temperature  replica  exchange  (TREX)  simulations  in  the  absence 
of  first  order  transitions/  '  “  the  situation  with  HREX  is  different. 

It  is  difficult  to  devise  a  universal  exchange  protocol  for  HREX 
simulations,  because  the  Hamiltonian  generally  depends  on  A 
nonlinearly.15,16  One  way  to  monitor  the  efficiency  of  HREX 
simulations  is  to  examine  the  number  of  round-trips  made  over 
the  simulation  time.42,71  A  replica  that  has  returned  to  its  initial  A 
after  visiting  both  A  =  0  and  A  =  1  states  in  either  order  is  said  to 
have  accomplished  a  round-trip.  Although  this  is  an  excellent 
measure  of  efficiency  in  theory,  simulations  might  not,  in 
practice,  be  long  enough  to  complete  even  a  single  round-trip. 
In  addition,  when  increasing  a  number  of  windows,  a  round-trip 
may  take  longer  time.  Therefore,  long  simulation  times  may  be 
required  to  use  this  measure  of  exchange  efficiency. 

Finally,  we  note  that  besides  the  A  =  1  window  in  SRTI,  scaling  of 
the  dihedral  and  soft-core  potentials  can  activate  hindered  transi¬ 
tions  in  several  neighboring  windows.  This  means  that  multiple 
windows  can  experience  enhanced  sampling.  While  a  standard 
round-trip  count  would  reflect  the  overall  efficiency  of  the  simula¬ 
tions,  it  may  not  capture  diffusion  of  the  conformations  from  all  of 
the  enhanced  windows  down  to  the  A  =  0  window.  Ultimately,  we 
gauge  the  sampling  gains  by  the  reduction  in  the  standard  deviations 
of  HREX-SRTI  compared  to  regular  SRTI,  which  is  an  independent 
measure  of  both  convergence  and  sampling  efficiency. 

3.32.2.  Acceptance  Ratio  Profile.  Alternatively,  one  can 
generate  an  acceptance  ratio  profile  for  each  pair  of  TI  windows 
that  is  adjacent  in  A  space  from  the  actual  simulation  data.15  A  flat 
profile  would  indicate  equally  probable  exchanges  between  adja¬ 
cent  windows  and,  consequently,  produce  the  largest  number  of 
round-trips  possible  for  any  given  time.  Therefore,  we  consider  the 
acceptance  ratio  profile  a  practical  alternative  to  the  round-trip 
count.  A  simple  inspection  of  the  acceptance  ratio  profile  could 
identify  problems  in  the  HREX  simulations.  The  pair  with  the 
lowest  ratio  across  all  of  the  window  pairs  limits  the  round-trip  rate. 

3.32.3.  Energy  Difference  Histograms.  The  acceptance  ratio 
profile  is  related  to  the  corresponding  double  and  single  energy 
difference  histograms.  Single  energy  difference  histograms  over  all 
adjacent  pairs  of  A,  and  AI+1  involve  the  corresponding  forward  and 
backward  energy  differences.  These  histograms  contain  valuable 
information  not  only  with  respect  to  the  efficiency  of  the  simula¬ 
tions  but  also  with  respect  to  their  validity.73  Furthermore,  the 
information  from  the  forward  and  backward  histograms  can  be 
used  to  estimate  the  free  energy  difference  between  the  adjacent 
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Figure  4.  Monitoring  efficiency  of  HREX-SRTI.  The  top  panel  illus¬ 
trates  the  acceptance  ratio  profile  for  indole  bound  to  the  L99A  mutant 
ofT4  lysozyme  with  the  V87p-A:Vlllp-A:L118p-A  reference  state.  The 
bottom  panel  shows  the  AA(/  histograms  in  color.  On  the  top  of  the 
histograms  is  the  actual  profile  of  the  mean  AAI;  value  (solid  black  line 
with  circles),  its  estimate  from  the  variance  (dashed  black  line  with 
diamonds),  and  the  contribution  of  the  derivative  of  the  <3V/9A> 
(dashed  black  fine  with  triangles),  which  are  derived  from  eq  4.  The 
HREX-SRTI  simulations  employed  23  windows,  each  run  for  2  ns  with 
p  =  2  soft-core  potentials  in  the  NPT  ensemble  at  1  atm  and  300  K. 


states,  though  such  an  estimator  may  be  suboptimal.43'73  In  cases 
with  linear  dependence  of  the  Hamiltonian  on  A,  the  single  energy 
difference  histogram  is  closely  related  to  the  corresponding 
histogram  of  the  energy  derivative  with  respect  to  A.43  The  double 
energy  difference  histogram  comprises  the  energy  change  of  the 
generalized  ensemble  that  enters  the  Metropolis  function  to 
decide  on  the  exchange  of  a  particular  pair.15  These  double  energy 
difference  histograms  are  perhaps  the  most  informative  for  the 
purpose  of  the  HREX.  They  could  be  computed  using  configura¬ 
tions  generated  by  regular  SRTI  for  a  given  set  of  A  values  at  an 
additional  expense  that  would  increase  the  cost  of  the  calculation 
to  that  of  HREX-SRTI. 

3.3.3.  Parameters  That  Influence  HREX  Efficiency.  Predicting 
the  dependence  of  the  double  energy  difference  or  the  accep¬ 
tance  ratio  profile  on  the  coupling  parameter  A  without  actually 
running  HREX  simulations  could  help  design  more  efficient 
simulations.  Indeed,  determining  an  optimal  set  of  A  values  that 
would  yield  a  uniform  acceptance  ratio  profile  will  maximize  the 
efficiency.  It  is  easy  to  show  that  the  mean  of  the  double  energy 
difference  histograms  and  the  (3V/3A)  are  related  according  to 
eq  4  (see  the  Appendix  for  derivation): 


(AA,j)  «/3(AA,j)2 


d  {  dV\ 
dA\3A/ 


(1,  +  A,-)/ 2y 


*/32(AA,;)2var(  ^ 


(A  +  X,)/2 


(4) 
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Figure  4  shows  an  overlay  of  the  A  A,,  distributions,  their  actual 
mean  value,  and  its  approximation  using  eq  4.  As  can  be  seen 
from  Figure  4,  eq  4  closely  approximates  the  mean  of  the  double 
energy  difference.  Furthermore,  the  contribution  of  the  mean 
second  derivative  is  significant  and  should  not  be  neglected. 
Thus,  eq  4  is  a  useful  starting  point  for  improving  HREX 
simulation  protocols. 

It  might  be  worthwhile  to  note  that  the  derivative  in  the  latter 
equation  is  related  to  the  variance  of  the  dV/dA:43'74 


(5) 


According  to  eq  5,  in  cases  when  the  potential  V  depends  on  A 
linearly,  the  variance  reduces  to  the  first  derivative  of  the  <3V/  9 A>. 
Because  the  variance  is  a  positive  quantity,  the  left-hand  side 
of  eq  4  would  always  be  positive.  Assuming  that  the  value  of  the 
mean  double  energy  difference  would  correspond  to  the  most 
probable  value  of  its  distribution,  one  might  expect  better 
acceptance  ratios  if  the  mean  is  closer  to  zero.  For  nonlinear 
dependence,  convexity  of  the  potential  with  respect  to  A  would 
play  an  important  role. 

3.3.3. 1.  Increasing  the  Number  of  Tl  Windows.  Equation  4 
suggests  that  reducing  the  spacing  between  adjacent  A’s 
(increasing  the  number  of  windows)  would  result  in  an  increase 
in  the  acceptance  probability.  Indeed,  similarly  to  our  previous 
study,22  increasing  the  number  of  windows  from  12  to  23  in 
the  system  with  the  V87-pA:Vlll-pA:L118-pA  reference 
state  considerably  improves  the  overall  acceptance  ratio  from 
less  than  10%  to  30%.  The  reduction  in  standard  deviation 
from  3.3  to  1.4  kcal/mol  attests  to  the  improved  efficiency. 
It  should  be  noted  that  the  cost  of  the  23-window  HREX-SRTI 
simulations  was  maintained  similar  to  that  of  the  12-window 
simulation  by  reducing  the  number  of  exchange  cycles  to 
1000. 

The  acceptance  ratio  profile  of  the  HREX-SRTI  simulations 
with  V87-pA:Vlll-pA:L118-pA  reference  and  23  windows 
(shown  in  Figure  4)  revealed  an  exchange  bottleneck.  Indeed, 
as  seen  in  Figure  4,  the  acceptance  ratio  for  the  pair  of  windows 
with  Aj  =  0.4091  and  Ai+1  =  0.4545  is  close  to  zero.  In  other  words, 
our  generalized  ensemble  was  divided,  with  replicas  sampling 
two  independent  regions  of  the  A  space.  This  prevented  the 
system  from  reaching  a  true  equilibrium.  Clearly,  the  23-window 
HREX  simulations,  while  greatly  improving  the  overall  accep¬ 
tance  ratio,  still  had  the  limitation  of  inefficient  exchanges  in  a 
specific  region  of  the  A  space. 

3. 3. 3. 2.  Position  of  Windows.  Although  generally  applicable, 
simply  increasing  the  number  of  windows  to  improve  the 
efficiency  of  HREX  simulations  is  not  a  viable  strategy. 
Certainly,  if  successive  TI  windows  were  kept  at  constant 
A  intervals,  this  strategy  would  not  remove  the  existing 
exchange  bottlenecks.  In  addition,  maintaining  the  number 
of  windows  at  the  same  level  as  in  regular  SRTI  simulations 
would  be  desirable  from  a  computational  cost  perspective. 
Therefore,  a  better  strategy  would  be  to  adjust  the  positions 
of  the  existing  TI  windows.  Importantly,  regular  TI  is 
sufficient  for  obtaining  the  dV/dA  variance  profile  that 
according  to  eq  4  could  be  used  to  determine  an  optimal 
set  of  A’s  for  efficient  HREX  simulations.  Alternatively, 
round-trip-based  methods  could  be  used  to  optimize  place¬ 
ment  of  A’s  as  has  been  done  in  the  context  of  simulated 
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Figure  5.  Optimization  of  HREX-SRTI  efficiency  by  altering  soft-core 
potential  parameters.  The  top  panel  illustrates  the  acceptance  ratio 
profiles  for  indole  bound  to  the  L99A  mutant  of  T4  lysozyme  with  the 
Vlllp-A  reference  state.  Unlike  in  Figure  4,  these  simulations  em¬ 
ployed  p  =  1  soft-core  potentials  with  three  different  values  for  the 
parameter  CL,  namely  0.3, 0.4,  and  0.5.  The  bottom  panel  shows  the  AAI; 
histograms  in  color  along  with  their  actual  mean  (solid  black  line  with 
circles),  the  estimate  of  the  mean  from  the  3V/3A  variance  (dashed 
black  line  with  diamonds),  and  the  contribution  of  the  derivative  of  the 
<3V/3A>  (dashed  black  line  with  triangles),  which  are  derived  from 
eq  4.  The  bottom  panel  only  shows  results  for  the  simulations  with  CL  = 
0.4  that  achieve  the  best  precision.  The  HREX-SRTI  simulations 
employed  only  12  windows,  each  run  for  4  ns  in  the  NPT  ensemble 
at  1  atm  and  300  K. 


scaling.42  However,  determining  the  optimal  set  of  A’s  is 
highly  system-dependent.  Therefore,  we  did  not  pursue  this 
strategy  here.  Fortunately,  eq  4  suggests  yet  another  strategy 
for  improving  HREX-SRTI  simulations. 

3.3.3.3.  Altering  Mean  Force  Profile.  Besides  manipulating 
the  number  and  position  of  TI  windows,  eq  4  indicates  that 
the  acceptance  probability  is  greatly  affected  by  the  shape  of  the 
dV/dA  variance  profile.  Therefore,  we  could  improve  the  ex¬ 
change  rates  by  changing  the  shape  of  the  profile  using  param¬ 
eters  we  have  at  our  disposal.  Specifically,  we  can  vary  parameters 
of  the  soft-core  potentials. 

Recall  the  general  form  of  the  soft-core  potentials  available  in 
GROMACS9,10,63'75  (eqs  1—3).  The  present  implementation  of 
GROMACS  supports  soft-cores  with  p  =  1  and  p  =  2.9,10  The  CL 
and  O  parameters  have  been  optimized  for  each  p  using  an  approach 
that  decouples  Coulomb  and  vdW  changes.76  For  the  SRTI  approach 
that  simultaneously  transforms  Coulomb  and  vdW  interactions,  these 
values  of  parameters  may  be  suboptimal. 

Soft-core  potentials  withp  =  1  are  well-suited  for  HREX-SRTI 
simulations.  Indeed,  our  earlier  SRTI  simulations  with  p  =  2 
potentials  had  the  limitation  of  sharp  peaks  in  the  (dV/  3 A)  profile 
near  A  =  0  as  previously  described  in  the  literature.9'67  These 
peaks  are  due  to  the  hydrogen  atoms  with  undefined  O  param¬ 
eters  and  present  an  obstacle  to  accurate  integration  within  the 
Fourier  beads  approach.22  In  order  to  integrate  these  sharp  peaks 
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Figure  6.  Comparison  of  the  HREX-SRTI  relative  binding  free  energy 
predictions  to  experimental  results.  The  binding  free  energies  are  relative 
to  benzene.  The  best-fitted  line  passing  through  the  origin  is  y  =  1.1  Ox, 
and  the  corresponding  correlation  coefficient  is  R~  =  0.76.  The  RESP 
HF/6-31G(d)  charge  model  was  used  for  the  cyclic  compounds,  and  the 
AM1BCC  model  was  used  for  the  acyclic  compounds.  With  the 
exception  of  indole,  all  of  the  simulations  employed  p  =  2  soft-cores 
with  default  parameters.  The  indole  prediction  incorporates  the  results 
obtained  with  the  optimized  p  =  1  soft-cores. 

properly,  additional  TI  windows  have  been  introduced.  In 
contrast,  the  profile  derived  using  p  =  1  soft-core  potentials  is 
devoid  of  such  peaks.10 

Switching  to  the  soft-core  potentials  with  p  =  1  improves  the 
acceptance  ratio  of  the  HREX-SRTI  simulations.  Indeed,  using 
default  parameters  with  only  12  windows,  we  obtained  modest 
acceptance  ratios  of  17%  in  the  system  that  had  the  V87-pA: 
VI 1  l-pA:Ll  18-pA  reference.  Recall  that  12-window  simulations 
with  p  =  2  failed  to  achieve  acceptance  ratios  above  10%. 
However,  as  discussed  above,  an  improved  acceptance  ratio  does 
not  guarantee  improved  efficiency.  The  analysis  of  the  accep¬ 
tance  ratio  profile  for  p  =  1  shows  that  windows  near  2  =  1 
exchange  poorly. 

3.33.4.  Optimization  of  Soft-Core  Parameters  for  HREX. 
The  optimization  of  the  CL  parameter  helps  to  flatten  the 
acceptance  ratio  profile,  rendering  the  HREX-SRTI  more  effi¬ 
cient.  Using  12-window  simulations  of  indole  with  the  V87-pA: 
VI  ll-pA:Ll  18-pA  reference,  we  varied  the  value  of  CL  from  its 
default  value  of  0.3  to  0.4  and  0.5  at  p  =  1.  Figure  5  summarizes 
the  acceptance  ratio  profiles  and  compares  predicted  and  actual 
mean  A  A;,  for  the  best  of  them.  While  at  a  =  0.5  we  obtained  the 
most  uniform  acceptance  ratio  profile,  the  overall  acceptance 
ratio  declined  to  9%  (Table  2).  At  the  intermediate  value  of  0.4, 
the  overall  acceptance  ratio  was  13%,  with  the  acceptance  ratio 
profile  demonstrating  sufficient  exchange  probabilities  near  A  =  1 . 
Clearly,  one  has  to  find  a  compromise  between  the  flatness  of  the 
acceptance  ratio  profile  and  the  overall  acceptance  ratio. 

To  choose  the  best  value  of  parameter  CL,  we  compared  the 
standard  deviations  from  each  set  of  simulations.  As  seen  in 
Table  2,  the  intermediate  value  CL  =  0.4  yielded  the  lowest 
standard  deviations  of  0.8  kcal/mol.  This  is  a  considerable 
improvement  for  the  system  with  the  triple  mutant  reference,  in 
comparison  to  the  regular  12-window  SRTI  simulation  with 
p  =  2,  which  had  a  standard  deviation  of  3.8  kcal/mol.  Inter¬ 
estingly,  simulations  with  CL  =  0.5  that  were  associated  with  an 
almost  flat  acceptance  ratio  profile  yielded  a  standard  deviation 
of  1.2  kcal/mol. 


3.33.5.  Transferability  of  Optimized  Parameters.  The  value 
of  CL  optimized  for  the  challenging  system  with  the  V87-pA: 
VI 1  l-pA:Ll  18-pA  reference  state  could  be  transferred  to 
improve  the  results  for  an  easier  system  with  the  Vlll-pA 
reference.  With  the  latter  reference  state,  only  the  VI 11  side 
chain  is  enhanced  along  with  the  ligand  in  HREX-SRTI.  Indeed, 
as  seen  in  Table  2,  for  p  =  2,  the  HREX-SRTI  simulations  with  12 
windows  and  the  Vlll-pA  reference  achieved  a  standard  devia¬ 
tion  of  1.8  kcal/mol.  In  fact,  HREX-SRTI  and  regular  SRTI  had 
identical  standard  deviations,  despite  the  fact  that  the  overall 
acceptance  ratio  of  the  former  was  19%.  With  23  windows,  both 
the  acceptance  ratio  and  the  standard  deviation  improved  to  48% 
and  0.5  kcal/mol,  respectively.  To  our  great  satisfaction,  the 
optimized  HREX-SRTI  with  p  =  1  demonstrated  significantly 
improved  convergence.  Remarkably,  we  were  able  to  achieve  an 
improved  standard  deviation  of  0.2  kcal/mol  using  only  12 
windows  with  an  overall  acceptance  ratio  of  29%. 

Even  with  the  optimized  protocol,  our  predictions  for  indole 
deviate  significantly  from  the  experimental  data.  Figure  6  shows 
the  comparison  of  our  predicted  binding  free  energies  with  the 
experimental  values  for  all  of  the  ligands,  including  the  optimized 
results  for  indole.  Indeed,  the  binding  free  energy  of  indole  using 
the  most  efficient  HREX-SRTI  protocol  is  still  overestimated  by 
1.5  kcal/mol. 

3.33.6.  Decoupling  of  Coulomb  and  vdW  Transforma¬ 
tions.  As  a  final  note  on  efficiency,  we  would  like  to  mention 
that  the  (dV/  92)  profile  would  change  radically  by  decoupling  of 
the  vdW  and  Coulomb  transformations.  Such  a  decoupling  is 
expected  to  render  the  overall  profile  smoother  and  hence 
improve  the  acceptance  ratio.21  Indeed,  HREX  simulations  with 
decoupled  vdW  and  Coulomb  transformations  were  recently 
reported.20,21  Unfortunately,  the  overall  acceptance  ratios  or 
acceptance  ratio  profiles  were  not  provided.  Since  the  HREX- 
SRTI  approach  was  originally  designed22  to  electrostatically 
guide  ligands  to  better  binding  poses,77’78  decoupling  of  the 
electrostatics  and  vdW  would  defeat  the  purpose. 

3.4.  HREX-SRTI  As  a  Conformational  Analysis  Tool.  It  is 
instructive  to  analyze  structural  transitions  in  the  L99A  T4 
lysozyme  mutant  in  complex  with  indole  in  an  attempt  to  explain 
the  observed  disagreement  with  experimental  values. 

3.4. 1.  Transitions  of  Hindered  Residues  Are  Indeed  Enhanced 
by  HREX-SRTI.  In  order  to  demonstrate  that  HREX-SRTI 
enhances  conformational  transitions  of  the  hindered  protein 
side  chains,  we  performed  additional  analyses  of  the  real  state 
trajectories  (2  =  0).  Figures  2  and  3  compare  the  dihedral  angles, 
described  in  the  Methods  section,  that  characterize  side  chain 
conformations  of  V87,  VI 11,  and  LI  18  where  applicable. 
Although  the  histograms  of  the  referred  to  dihedral  angles  may 
look  similar,  the  time  series  clearly  demonstrate  the  increased 
number  of  transitions  in  HREX-SRTI.  In  most  cases,  each  side 
chain  samples  multiple  conformational  basins,  supporting  our 
earlier  conclusion  regarding  their  contributions  to  the  binding 
free  energies.  Interestingly,  the  V87  side  chain  appears  to  con¬ 
sistently  prefer  a  particular  conformation. 

3.4.2.  Indole  Transitions  in  the  Binding  Pocket.  We  also 
verified  that  indole  could  flip  its  benzene  plane  in  HREX-SRTI 
simulations  with  triple  mutant  reference  states.  The  analysis  of 
the  real  state  trajectories  revealed  multiple  indole  orientations  in 
the  binding  pocket.  In  one  of  the  orientations,  the  NH  group  of 
indole  persistently  hydrogen-bonded  with  the  S  atom  of  M102. 
In  another  orientation,  the  same  NH  group  was  involved  in 
hydrogen  bonding  with  the  backbone  carbonyl  group  of  V87. 
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In  the  HREX-SRTI  simulations  with  the  single  mutant  VI 1  lp-A 
reference  state,  the  ligand  in-plane  rotation  was  enhanced,  but 
the  plane  flipping  transitions  remained  hindered. 

3.4.3.  Water  Does  Not  Bind  T4  Lysozyme  Mutant  with  Indole. 
The  indole  disagreement  with  experimental  values  is  not  attributable 
to  a  lack  of  water  access  to  the  binding  site.  Although  our  simulations 
with  the  VI 1  lp-A  reference  state  are  not  designed  to  open  the 
normally  sealed  active  site,  we  have  tested  this  hypothesis  using  the 
V87p-A:Vlllp-A:L118p-A  reference  state.  By  inspecting  the  MD 
trajectories,  we  observed  that  water  molecules  did  penetrate  the 
pocket  of  the  triple  mutant  reference  state  (A  =  l).  However,  none 
of  the  configurations  with  water  molecules  inside  the  binding  site 
reached  the  real  state  (A  =  0)  during  the  HREX-SRTI  simulations. 
This  strongly  suggests  that  water  is  not  responsible  for  the  observed 
discrepancy. 


The  generalized  ensemble  average  value  of  the  double  differ¬ 
ence  for  a  given  pair  of  adjacent  replicas  is  then 

<AA,j)  =  /3(AVtj)^  +  /3<AV;,)a.  (A3) 


Using  second  order  Taylor  expansion: 
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where 


A  A^  —  Xj  —  A,-  (A5) 

we  obtain  the  following  expression,  which  is  equivalent  to  eq  4  in 
the  main  text: 


4.  CONCLUSIONS 

This  study  presented  a  practical  application  of  the  SRTI 
approach  to  compute  relative  binding  free  energies  of  small 
molecules  to  a  challenging  binding  site  in  the  L99A  mutant 
of  T4  lysozyme.  With  the  HREX  option,  SRTI  successfully 
enhanced  sampling  of  the  hindered  transitions  of  protein  side 
chains  and  bound  ligands.  Achieving  efficient  HREX  simulations 
improves  the  quality  of  predictions.  However,  the  commonly 
used  overall  acceptance  ratio  is  not  a  good  indicator  of  the 
efficiency  of  the  HREX.  Instead,  acceptance  ratio  profiles  should 
be  examined  and  whenever  possible  made  uniform  by  adjusting 
simulation  parameters.  To  aid  the  future  design  of  efficient 
simulation  protocols,  we  have  provided  a  useful  relationship 
between  the  mean  exchange  energy  and  the  corresponding  3V/ 
3A  variance  profiles.  Guided  by  the  relationship,  we  demon¬ 
strated  that  judicial  changes  in  the  soft-core  potentials  consider¬ 
ably  improved  HREX-SRTI  simulation  efficiency.  Overall,  the 
HREX-SRTI  predicted  relative  binding  free  energies  for  a  series 
of  12  ligands  with  an  RMSE  of  0.3  kcal/mol  comparable  to 
experimental  data.  Ultimately,  improving  efficiency  of  the  HREX 
simulations  may  further  reduce  computational  cost  and  increase 
the  precision  of  the  predictions. 

Note,  while  this  paper  was  under  revision,  we  discovered  a 
paper  by  Steiner  and  coauthors  that  used  an  approach  identical  to 
HREX-SRTI  to  compute  relative  free  energies  of  a  number  of 
ligands  to  Plasmepsin  II,  albeit  without  accelerating  any  residues 
of  the  protein.7'1 

■  APPENDIX 

For  a  configuration  R,  the  vertical  excitation  energy 
Aj;(R)  from  a  Hamiltonian  at  A,-  to  that  at  A;  is  calculated  as 
follows: 


A V^R)  =  V(R,Xj)  -  V (R,  A,)  (Al) 

The  overall  potential  energy  change  AAij(R,Rl)  for  the 
generalized  ensemble  upon  Hamiltonian  exchange  between 
two  configurations  R  and  R'  from  the  adjacent  windows  A, 
and  kj  is 

AA;;-(R,R')  =  PlAViMx,  +  /3[AV),(R')]1;  (A2) 

The  subscripts  after  the  square  brackets  indicate  the  Hamiltonian 
of  the  simulations  used  to  obtain  the  respective  configurations. 
Only  two  configurations  representing  each  A  are  involved. 


<AA1?>  « /3  AX, 


3A2  /  A  dA  \  3A  i  , 


(«+/)/V 


(A6) 


Recalling  eq  5  from  the  main  text  we  obtain  the  final  relation. 

q2  /  a  1  \2  /SV\ 

~rctr(  -  ) 

3A  /; 

'  V>)/2 


(AAjj)  »  /32(AA,j)2  var 


(A7) 


■  AUTHOR  INFORMATION 

Corresponding  Author 

*E-mail :  ikhavrutskii  (3)bioanalysis .  org. 


■  ACKNOWLEDGMENT 

We  would  like  to  thank  Dr.  In-Chul  Yeh,  Dr.  Michael  S.  Lee, 
and  Dr.  Hyung-June  Woo  for  stimulating  scientific  discussions. 
Also,  we  acknowledge  the  National  Cancer  Institute  (NCI)  for 
an  allocation  of  computing  time  and  staff  support  at  the 
Advanced  Biomedical  Computing  Center  (ABCC)  at  National 
Cancer  Institute,  Frederick,  Maryland.  This  work  was  sponsored 
by  the  U.S.  Department  of  Defense  High  Performance  Comput¬ 
ing  Modernization  Program  (HPCMP)  under  the  High  Perfor¬ 
mance  Computing  Software  Applications  Institutes  (HSAI) 
initiative.  Disclaimer:  The  opinions  and  assertions  contained 
herein  are  the  private  views  of  the  authors  and  are  not  to  be 
construed  as  official  or  as  reflecting  the  views  of  the  U.S.  Army  or 
the  U.S.  Department  of  Defense. 

■  REFERENCES 

(1)  Boyce,  S.  E.;  Mobley,  D.  L.;  Rocklin,  G.  J.;  Graves,  A.  P.;  Dill, 
K.  A.  /.  Mol.  Biol.  2009,  394,  747. 

(2)  Gallicchio,  E.;  Levy,  R  M.  Curr.  Opin.  Struct.  Biol.  2011, 21,  161. 

(3)  Aleksandrov,  A.;  Thompson,  D.;  Simonson,  T.  /.  Mol.  Recognit. 
2010,  23,  117. 

(4)  Bruckner,  S.;  Boresch,  S. /.  Comput.  Chem.  2011,  32,  1303. 

(5)  Chodera,  J.  D.;  Mobley,  D.  L.;  Shirts,  M.  R;  Dixon,  R  W.; 
Branson,  K,;  Pande,  V.  S.  Curr.  Opin.  Struct.  Biol.  2011,  21,  150. 

(6)  Leavitt,  S.;  Freire,  E.  Curr.  Opin.  Struct.  Biol.  2001,  11,  560. 

(7)  Chaires,  J.  B.  Annu.  Rev.  Biophys.  2008,  37,  135. 

(8)  Straatsma,  T.  P.;  McCammon,  J.  A  ].  Chem.  Phys.  1994, 101,  5032. 


3010 


dx.doi.org/10.1021/ct2003786  | J.  Chem.  Theory  Comput.  201 1,  7,  3001-301 1 


Journal  of  Chemical  Theory  and  Computation 


ARTICLE 


(9)  Beutler,  T.  C.;  Mark,  A.  E.;  van  Schaik,  R  C.;  Gerber,  P.  R;  van 
Gunsteren,  W.  F.  Chem.  Phys.  Lett.  1994,  222,  529. 

(10)  Shirts,  M.  R;  Pande,  V.  S. /.  Chem.  Phys.  2005,  122,  134508. 

(11)  Straatsma,  T.  P.;  McCammon, J.  A./.  Chem.  Phys.  1991, 95, 1 175. 

(12)  Kwak,  W.;  Hansmann,  U.  H.  E.  Phys.  Rev.  Lett.  2005,  95, 138102. 

(13)  Fukunishi,  H.;  Watanabe,  O.;  Takada,  S.  /.  Chem.  Phys.  2002, 
116,  9058. 

( 14)  Sugita,  Y.;  Kitao,  A.;  Okamoto,  Y.  /.  Chem.  Phys.  2000, 113, 6042. 

(15)  Hritz,  J.;  Oostenbrink,  C .  J.  Chem.  Phys.  2007,  127,  204104. 

(16)  Hritz,  J.;  Oostenbrink,  C .  J.  Chem.  Phys.  2008,  128,  144121. 

(17)  Hritz,  J.;  Oostenbrink,  C. /.  Phys.  Chem.  B  2009,  113,  12711. 

(18)  Woods,  C.  J.;  Essex,  J.  W.;  King,  M.  A  /.  Phys.  Chem.  B  2003, 
107,  13703. 

(19)  Woods,  C.  J.;  King,  M.  A;  Essex,  J.  W.  Lect.  Notes  Comput.  Sci. 
Eng.  2006,  49,  251. 

(20)  Jiang,  W.;  Hodoscek,  M.;  Roux,  B.  J.  Chem.  Theory  Comput. 
2009,  5,  2583. 

(21)  Jiang,  W.;  Roux,  B. /.  Chem.  Theory  Comput.  2010,  6,  2559. 

(22)  Khavrutskii,  I.  V.;  Wallqvist,  A.  /.  Chem.  Theory  Comput.  2010, 
6,  3427. 

(23)  Shirts,  M.  R;  Pande,  V.  S.  /.  Chem.  Phys.  2005,  122,  144107. 

(24)  Shirts,  M.  R;  Chodera,J.  D ./.  Chem.  Phys.  2008,  129,  124105. 

(25)  Shirts,  M.  R;  Bair,  E.;  Hooker,  G.;  Pande,  V.  S.  Phys.  Rev.  Lett. 

2003,  91,  140601. 

(26)  Bennett,  C.  H.  /.  Comput.  Phys.  1976,  22,  245. 

(27)  Kumar,  S.;  Bouzida,  D.;  Swendsen,  R  H.;  Kollman,  P.  A.; 
Rosenberg,  J.  M. /.  Comput.  Chem.  1992,  13,  1011. 

(28)  Deng,  Y.;  Roux,  B.  J.  Chem.  Theory  Comput.  2006,  2,  1255. 

(29)  Mobley,  D.  L.;  Chodera,  J.  D.;  Dill,  K.  A.  J.  Chem.  Theory 
Comput.  2007,  3,  1231. 

(30)  Mobley,  D.  L.;  Graves,  A.  P.;  Chodera,  J.  D.;  McReynolds,  A.  C.; 
Shoichet,  B.  K.;  Dill,  K.  A  /.  Mol.  Biol.  2007,  371,  1118. 

(31)  Mark,  A.  E.;  van  Gunsteren,  W.  F.;  Berendsen,  H.  J.  C.  J.  Chem. 
Phys.  1991,  94,  3808. 

(32)  Wereszczynski,  J.;  McCammon,  J.  A.  /.  Chem.  Theory  Comput. 
2011,  6,  3285. 

(33)  Fajer,  M.;  Hamelberg,  D.;  McCammon,  J.  A.  /.  Chem.  Theory 
Comput.  2008,  4,  1565. 

(34)  Fajer,  M.;  Swift,  R  V.;  McCammon,  J.  A.  J.  Comput.  Chem. 
2009,  30,  1719. 

(35)  Hamelberg,  D.;  Mongan,  J.;  McCammon,  J.  A.  /.  Chem.  Phys. 

2004,  120,  11919. 

(36)  Kong,  X.  J.;  Brooks,  C.  L.  J.  Chem.  Phys.  1996,  105,  2414. 

(37)  Banba,  S.;  Guo,  Z.;  Brooks,  C.  L ./.  Phys.  Chem.  B  2000, 104,  6903. 

(38)  Bitetti-Putzer,  R;  Yang,  W.;  Karplus,  M.  Chem.  Phys.  Lett.  2003, 
377,  633. 

(39)  Eriksson,  M.  A.  L.;  Pitera,  J.;  Kollman,  P.  A  /.  Med.  Chem.  1999, 
42,  868. 

(40)  Pitera,  J.;  Kollman,  P .  J.  Am.  Chem.  Soc.  1998,  120,  7557. 

(41)  Li,  H.;  Fajer,  M.;  Yang,  W.  /.  Chem.  Phys.  2007,  126,  024106. 

(42)  Zheng,  L.;  Yang,  W.  J.  Chem.  Phys.  2008,  129,  124107. 

(43)  Min,  D.;  Li,  H.;  Li,  G.;  Bitetti-Putzer,  R;  Yang,  W ./.  Chem.  Phys. 
2007,  126,  144109. 

(44)  Sadowski,  J.;  Gasteiger,  J.;  Klebe,  G.  /.  Chem.  Inf.  Comput.  Sci. 
1994,  34,  1000. 

(45)  Renner,  S.;  Schwab,  C.  H.;  Gasteiger,  J.;  Schneider,  G. /.  Chem. 
Inf.  Model.  2006,  46,  2324. 

(46)  Wang,  J.;  Wolf,  R  M.;  Caldwell,  J.  W.;  Kollman,  P.  A;  Case, 
D.  A.  /.  Comput.  Chem.  2004,  25,  1157. 

(47)  Dewar,  M.  J.  S.;  Zoebisch,  E.  G.;  Healy,  E.  F.;  Stewart,  J.  J.  P. 
J.  Am.  Chem.  Soc.  1985,  107,  3902. 

(  48  )  Stewart,  J.  J.  P.  MOPAC7 ;  University  of  T exas,  Austin :  Austin,  TX 

(49)  Jakalian,  A;  Bush,  B.  L.;  Jack,  D.  B.;  Bayly,  C.  I./.  Comput.  Chem. 
2000,  21,  132. 

(50)  Jakalian,  A.;  Jack,  D.  B.;  Bayly,  C.  I.  /.  Comput.  Chem.  2002, 
23,  1623. 

(51)  Wang,  J.  Antechamber  2009,  1,  2. 


(52)  Wang,  J.;  Wang,  W.;  Kollman,  P.  A;  Case,  D.  A.  J.  Mol.  Graphics 
Model.  2006,  25,  247. 

(53)  Nicholls,  A.;  Mobley,  D.  L.;  Guthrie,  J.  P.;  Chodera,  J.  D.;  Bayly, 
C.  I.;  Cooper,  M.  D.;  Pande,  V.  S.  /.  Med.  Chem.  2008,  51,  769. 

(54)  Mobley,  D.  L.;  Bayly,  C.  I.;  Cooper,  M.  D.;  Dill,  K.  A.  J.  Phys. 
Chem.  B  2009,  113,  4533. 

(55)  Mobley,  D.  L.;  Bayly,  C.  I.;  Cooper,  M.  D.;  Shirts,  M.  R;  Dill, 
K.  A.  /.  Chem.  Theory.  Comput.  2009,  5,  350. 

(56)  Shivakumar,  D.;  Deng,  Y.j  Roux,  B ./.  Chem.  Theory  Comput.  2009, 
5,  919. 

(57)  Frisch,  M  J.;  Trucks,  G.  W.;  Schlegel,  H.  B.;  Scuseria,  G.  E.;  Robb, 
M.  A;  Cheeseman,  J.  R;  Scalmani,  G.;  Barone,  V.;  Mennucci,  B.;  Petersson, 

G.  A;  Nakatsuji,  H.;  Caricato,  M.;  Li,  X;  Hratchian,  H.  P.;  Izmaylov,  A  F.; 
Bloino,  J.;  Zheng,  G.;  Sonnenberg,  J.  L.;  Hada,  M.;  Ehara,  M.;  Toyota,  K; 
Fukuda,  R;  Hasegawa,  J.;  Ishida,  M;  Nakajima,  T.;  Honda,  Y.;  Kitao,  O.;  Nakai, 

H. ;  Vreven,  T.j  Montgomery,  J.  A,  Jr.;  Peralta,  J.  E.;  Ogliaro,  F.;  Bearpark,  M; 
Heyd,  J.  J.;  Brothers,  E.;  Kudin,  K  N.;  Staroverov,  V.  N.;  Kobayashi,  R; 
Normand,J.;  Raghavachari,  K;  Rendell,  A;  Burant,J.  C.;  Iyengar,  S.  S.;  Tomasi, 
J.;  Cossi,  M;  Rega,  N.;  Millam,  N.  J.;  Klene,  M.;  Knox,  J.  E.;  Cross,  J.  B.;  Bakken, 
V.;  Adamo,  C.;  Jaramillo,  J.;  Gomperts,  R;  Stratmann,  R  E.;  Yazyev,  O.;  Austin, 
A  J.;  Cammi,  R;  Pomelli,  C.;  Ochterski,  J.  W.;  Martin,  R  L.;  Morokuma,  K; 
Zakrzewski,  V.  G.;  Voth,  G.  A;  Salvador,  P.;  Dannenberg,  J.  J.;  Dapprich,  S.; 
Daniels,  A  D.;  Farkas,  O.;  Foresman,  J.  B.;  Ortiz,  J.  V.;  Cioslowski,  J.;  Fox,  D.  J. 
Gaussian  09,  Revision  Al;  Wallingford,  CT,  2009. 

(58)  Bayly,  C.  I.;  Cieplak,  P.;  Cornell,  W.;  Kollman,  P.  A.  J.  Phys. 
Chem.  1993,  97,  10269. 

(59)  Morton,  A.;  Matthews,  B.  W.  Biochemistry  1995,  34,  8576. 

(60)  Bron,  C.;  Kerbosch,  J.  Commun.  ACM  1973,  16,  575. 

(61)  Hornak,  V.;  Abel,  R;  Okur,  A.;  Strockbine,  B.;  Roitberg,  A.; 
Simmerling,  C.  Proteins  2006,  65,  712. 

(62)  Mobley,  D.  L.;  Chodera,  J.  D.;  Dill,  K.  A.  /.  Chem.  Phys.  2006, 
125,  084902. 

(63)  Van  Der  Spoel,  D.;  Lindahl,  E.;  Hess,  B.;  Groenhof,  G.;  Mark, 
A  E.;  Berendsen,  H.  J.  /.  Comput.  Chem.  2005,  26,  1701. 

(64)  Berendsen,  H.  J.  C.;  van  der  Spoel,  D.;  van,  D.  R  Comput.  Phys. 
Commun.  1995,  91,  43. 

(65)  Hess,  B.;  Kutzner,  C.;  van  der  Spoel,  D.;  Lindahl,  E.  J.  Chem. 
Theory  Comput.  2008,  4,  435. 

(66)  Van  Der  Spoel,  D.;  Lindahl,  E.;  Hess,  B.;  Groenhof,  G.;  Mark, 
A.  E.;  Berendsen,  H.  J.  C .  J.  Comput.  Chem.  2005,  26,  1701. 

(67)  van  der  Spoel,  D.;  Lindahl,  E.;  Hess,  B.;  Kutzner,  C.;  van 
Buuren,  A  R;  Apol,  E.;  Meulenhoff,  P.  J.;  Tieleman,  D.  P.;  Sijbers,  A  L.  T. 
M.;  Feenstra,  R  A;  Drunen,  R  v.;  Berendsen,  H.  J.  C.  GROMACS  User 
Manual  Version  4.0-,  The  GROMACS  development  team:  Groningen, 
The  Netherlands,  2006. 

(68)  Blondel,  A.  J.  Comput.  Chem.  2004,  25,  985. 

(69)  Hess,  B.;  Bekker,  H.;  Berendsen,  H.  J.  C.;  Fraaije,  J.  G.  E.  M. 
/.  Comput.  Chem.  1997,  18,  1463. 

(70)  Khavrutskii,  I.  V.;  Fajer,  M.;  McCammon,  J.  A.  J.  Chem.  Theory 
Comput.  2008,  4,  1541. 

(71)  Trebst,  S.;  Troyer,  M.;  Hansmann,  U.  H.  E. /.  Chem.  Phys.  2006, 
124,  174903. 

(72)  Yeh,  I.-C.;  Lee,  M.  S.;  Olson,  M.  A  /.  Phys.  Chem.  B  2008, 
112,  15064. 

(73)  Pohorille,  A;  Jarzynski,  C.;  Chipot,  C.  /.  Phys.  Chem.  B  2010, 
114,  10235. 

(74)  Simonson,  T.  Mol.  Phys.  1993,  80,  441. 

(75)  Lange,  O.  F.;  Schaefer,  L.  V.;  Grubmueller,  H./.  Comput.  Chem. 
2006,  27,  1693. 

(76)  Bash,  P.  A.;  Singh,  U.  C.;  Langridge,  R;  Kollman,  P.  A.  Science 
1987,  236,  564. 

(77)  Sines,  J.  J.;  Allison,  S.  A;  McCammon,  J.  A  Biochemistry  1990, 
29,9403. 

(78)  Tan,  R  C.;  Truong,  T.  N.;  McCammon,  J.  A;  Sussman,  J.  L. 
Biochemistry  1993,  32,  401. 

(79)  Steiner,  D.;  Oostenbrink,  C.;  Diederich,  F.;  Zurcher,  M.;  van 
Gunsteren,  W.  F.  /.  Comput.  Chem.  2011,  32,  1801. 


3011 


dx.doi.org/10.1021/ct2003786  | J.  Chem.  Theory  Comput.  201 1,  7,  3001-301 1 


